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We develop a time-dependent variational Monte Carlo (t-VMC) method for quantum dynamics 
of strongly correlated electrons. The t-VMC method has been recently applied to bosonic systems 
and quantum spin systems. Here, we propose a time-dependent trial wave function with many 
variational parameters, which is suitable for nonequilibrium strongly correlated electron systems. 
As the trial state, we adopt the generalized pair-product wave function with correlation factors 
and quantum-number projections. This trial wave function has been proven to accurately describe 
ground states of strongly correlated electron systems. To show the accuracy and efficiency of our 
trial wave function in nonequilibrium states as well, we present our benchmark results for relaxation 
dynamics during and after interaction quench protocols of fermionic Hubbard models. We find that 
our trial wave function well reproduces the exact results for the time evolution of physical quantities 
such as energy, momentum distribution, spin structure factor, and superconducting correlations. 
These results show that the t-VMC with our trial wave function offers an efficient and accurate way 
to study challenging problems of nonequilibrium dynamics in strongly correlated electron systems. 


I. INTRODUCTION 

Quantum systems with strong many-body correla¬ 
tions in equilibrium show intriguing properties such 
as the metal-insulator transition^ and high-temperature 
superconductivity^i^. Recently, because of potential 
routes to realizing intriguing phenomena that are not at¬ 
tainable in the equilibrium, strongly correlated electron 
systems driven out of equilibrium have attracted much 
attention. In fact, owing to the development of experi¬ 
mental techniques, we have been able to control or realize 
unprecedented phases and their phase transitions by ap¬ 
plying strong and short pulse of external fields such as 
intensive laser pumping^^— . 

For satisfactory theoretical understanding of quantum 
dynamics of many-body systems, we have to solve many- 
body time-dependent Schrddinger equation = 

'H{t) \tp{t)), where and 'H(t) represent the wave 

function and the Hamiltonian at time t, respectively. 
The formal solution of the time-dependent Schrddinger 
equation is given by |V'(0) = Texp(—f 'H{s)ds) |'!/'(0))- 
Here, T represents the time ordering. However, such an 
approach is tractable only for small many-body systems, 
because the Hilbert space grows exponentially as the sys¬ 
tem size increases. Furthermore, reduction to an effec¬ 
tive single-particle problems such as the time-dependent 
Hartree-Fock methodic does not give us accurate re¬ 
sults. To treat larger systems accurately, there exist sev¬ 
eral numerical methods such as time-dependent density 
matrix renormalization group method(DMRG)^“— , and 
nonequilibrium dynamical mean field theory (DMFT)I^ 
and quantum Monte Carlo(QMC) methodic. However, 
DMRG and DMFT have difficulties in treating large sys¬ 
tems in two or three spatial dimensions when one wishes 
to treat spatial correlations and fluctuations accurately. 
In order to include non-local correlations, the dynami¬ 
cal cluster approximation(DGA) has been proposed as 


an extension of the DMFT—. However, it requires high 
computational costs for a large cluster size. Although the 
QMG method can treat hnite systems exactly, its applica¬ 
tions are very limited due to the notorious negative-sign 
problem. Recently, to overcome the above difficulties, 
Garleo et al. developed the time-dependent variational 
Monte Garlo (t-VMG) method and optimized Jastrow 
factors in bosonic systems2ii^. This method can be for¬ 
mulated based on the time dependent variational prin¬ 
ciple (TDVP)^!^^!^!. The t-VMG method has been ap¬ 
plied not only to bosonic systems2ii2^ but also to spin 
models^. However, to the best of our knowledge, its 
application to correlated electron systems has not been 
successful yet. This limitation may be ascribed to the 
difficulty of constructing an accurate trial wave function 
for such systems. Therefore, the proposal of accurate 
trial wave functions suitable for nonequilibrium electron 
systems in the t-VMG method is desirable. 

The purpose of this study is to propose an accurate 
and efficient trial wave function for strongly correlated 
electron systems out of equilibrium in the t-VMC frame¬ 
work. To achieve this purpose, we focus on highly ac¬ 
curate trial wave functions for ground states of corre¬ 
lated electron systems. In such systems, many studies 
have attempted to construct an accurate and efficient 
trial wave function in the variational Monte Carlo (VMC) 
framework^^^— . In Ref.[3l|, Tahara and one of the au¬ 
thors have reduced the biases by using the quantum- 
number projections and introducing many variational pa¬ 
rameters to one-body part. They have adopted a gener¬ 
alized pair-product wave function as a one-body part be¬ 
cause it can flexibly describe different competing phases 
such as correlated metals, antiferromagnetic states and 
superconducting states. This improved trial wave func¬ 
tion has proven to be highly accurate for ground states of 
strongly correlated electron systems^i^i^. In this paper, 
we show that this trial wave function is an accurate and 
efficient one even for nonequilibrium strongly correlated 
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electron systems in the t-VMC framework. 

The organization of this paper is as follows. In Sec. I 
I, we introduce the TDVP which enables us to obtain an 
optimal time-dependent trial wave function and formu¬ 
late the t-VMC method by using the TDVP. Section II 
I describes a trial wave function with a large number of 
variational parameters for nonequilibrium strongly corre¬ 
lated electron systems. In Sec. IV, we show the accuracy 
and efficiency of our trial wave function by presenting 
several benchmark results. Finally, we summarize our 
work in Sec. V. 


II. TIME-DEPENDENT VARIATIONAL 
PRINCIPLE 


evolution of the variational parameters^ii^SiiHi^: 

N„ 


Oik = 


dok 

dt 


= %igu 


(3) 


where a matrix S and a vector g are described as 

Ski = (OlOi) - (Ol) (Oi) , (4) 

gk = (Ol-H) - (Ol) in) , (5) 

respectively. In the t-VMC method, we estimate an ex¬ 
pectation value (A) = by the Markov-chain 

Monte Carlo method. The derivative operators Ok and 
ol are defined by using real space configurations of elec¬ 
trons {a:} as 


The time-dependent variational principle (TDVP) pro¬ 
posed by McLachlan is a variational principle for time- 
dependent wave functions2^. In this principle, we con¬ 
sider a distance between i-^ \4’ct) and T-L \ tpa) where a = 
{afc|/c = I, • • • , Np} represent time-dependent variational 
parameters. By definition, the distance satisfies the in¬ 
equality 


min 

OL 


\tpa) -'H\lpa) 


> 0 , 


( 1 ) 


O/c = ^ lx) Ok{x) (x|, 

=^l^) Ol{x) (x|. 


respectively. Here, 
Ok{x) 
Ol{x) 


1 d 

(x|V^a) da 
1 d 

{iI)ol\x) da 


k 


* 

k 


(V'ala;)- 


( 6 ) 


(7) 


where the equality holds if \4’a) is the solution of the 
time-dependent Schrodinger equation. Here, the norm 
II 14') II is defined as the square root of an inner product of 
a wave function |'I'), i.e. || Idt) || = ^ ('I'I'I'). If we could 
optimize the variational parameters at each time step 
such that the equality holds, we obtain the exact solution 
of IV’ct). If a trial wave function well approximates the 
exact solution of the time-dependent Schrodinger equa¬ 
tion, the value of the lower bound should be small. Based 
on this idea, we optimize variational parameters at each 
time-step such that the distance is minimized. Origi¬ 
nally, the TDVP was applied in the field of quantum 
chemistry24i2^. Recently, a similar principle has been 
applied to the matrix product state for quantum spin 
models^i^, the bosonic Jastrow-type wave function for 
the Bose-Hubbard mode P^'^^ and the Gutzwiller approx¬ 
imation for strongly correlated electron systems^^— . 

Although exact time evolution is unitary, and thus, the 
norm {ipa\'4’a) is conserved, it is not necessary conserved 
in TDVP [Eq. ((T|)]. To remove the restriction on the 
norm, we use a TDVP for norm-independent dynamics^^. 


The differential equation [Eq. ([3])] is called TDVP 
equatior^^^. This TDVP equation can also be derived 
by minimizing the time-dependent actio n^^'^° . If we use 
the time-dependent variational principle for imaginary 
time evolution t = —ir and solve TDVP equation by us¬ 
ing Euler method, we obtain the stochastic reconfigura¬ 
tion scheme proposed by Sollera3^. The TDVP equation 
has a symplectic property^iiiSiil. This property leads 
to the energy conservation if the Hamiltonian is time- 
independent and we could calculate the derivative of pa¬ 
rameters dk exactly. 

In this study, in order to solve the TDVP equation, 
we use the fourth-order Runge-Kutta method which pro¬ 
vides us with a stable and efficient way to perform the 
time integration. Note that the Runge-Kutta method is 
not a symplectic integral method. Eurthermore, there 
are stochastic errors in the Monte Carlo calculation of 
quantities such as {OItL) and {OlOi). These cause the 
breaking of the symplectic property of the TDVP equa¬ 
tion. Nevertheless, we observed that the energy is con¬ 
served with high accuracy as we show in Sec. IV. 


min 

q: 


V (V'al^/’a) ) 


Hoc) -'HHoc) 


> 0 .( 2 ) 


III. TIME-DEPENDENT VARIATIONAL WAVE 
FUNCTION 


The details of the TDVP for norm-independent dynam¬ 
ics is described in Appendix [^ Based on this TDVP, we 
can derive the differential equation of the time-dependent 
variational parameters. Namely, by solving the mini¬ 
mization problem on the distance we obtain the time 


In the t-VMC method, the choice of trial wave func¬ 
tions is important. As a trial wave function, we adopt 
the form of 


Hit)) =cp{t)\(i){t)), 


( 8 ) 
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which has been used for equilibrium systems in Refs. mi¬ 
ls^. Here, £ represents quantum-number projections 
which recovers the symmetries the wave function should 
have throughout the time evolution, and V{t) represents 
correlation factors. For the one-body part \(j){t)), we em¬ 
ploy the pair-product wave function. In addition, we 
include backflow correlations in the pair-product wave 
function for lattice mode P^d^ . The time-dependent vari¬ 
ational parameters are included in the correlation factors 
as well as in the one-body part with the backflow corre¬ 
lations. Note that these variational parameters should 
be treated as complex numbers because the variational 
parameters evolve as complex numbers in the present 
method. In this section, we describe each component 
in detail. 


A. One-body Part 

In the conventional VMC, the Slater determinant with 
small variational parameters is used as the one-body part. 
In order to improve the conventional one-body part, we 
assume the form of the pair-product wave function with 
many variational parameters^: 

N/2 

|0), (9) 

where N is the number of electrons, Ng is the system 
size, (cicr) is a creation (annihilation) operator of an 
electron with spin a on the site £ and the pairing ampli¬ 
tude fij is treated as variational parameters. This pair- 
product wave function is a general form of a Hartree- 
Fock-Bogoliubov-type wave function which allows anti¬ 
ferromagnetic and superconducting order o^^d^d^^ . Thus, 
it takes an advantage of flexibly describing the paramag¬ 
netic metals, antiferromagnetic ordered states and super¬ 
conducting states with any type of frequency independent 
gap on equal footing. 


/ N, 

!</>)= /bcItSd 


In the limit g —>■ oo, the Gutzwiller wave function Vg {(f) 
corresponds to a state which contains no double occupa¬ 
tion. The Gutzwiller factor is a simple way to improve 
a mean-field wave function such as a Slater determinant 
and the pair-product wave function. However, it was nu¬ 
merically proven that the Gutzwiller wave function can¬ 
not describe the nonmagnetic Mott transition in any fi¬ 
nite dimensional systemsi^. The main reason for this 
is that the Gutzwiller factor only includes the on-site 
correlation. Although some doubly-occupied (doublon) 
and empty (holon) sites exit in the Mott insulator where 
charge fluctuations are allowed, the doublon and holon 
have to be bound in realizing an insulating behavior. 

In order to describe the Mott transition, the Jastrow 
factor is introduced^!: 

Pj = exp - l){nj - I)j , (II) 

where n* = rii^ -|- 714 , U 4 = v{ri - rj) = v{rj - n) 
are the variational parameters, and ri[rj) represents the 
position vector of the site i{j). The Jastrow factor can be 
represented by using a doublon number operator and a 
holon number operator. The electron number operator ni 
is written as rii = I + Di — Hi, where Di = 714714 is the 
doublon number operator and Hi = {1 — 774 ) (1 — 774 ) 
is the holon number operator. By using this relation, 
Jastrow factor becomes 

Pj = exp ^ Vij{DiDj -I- HiHj - DiHj - HiDj^ . 

( 12 ) 

Equation m shows that the Jastrow factor includes 
off-site repulsive correlations between doublon-doublon 
(holon-holon) pairs and attractive correlations between 
doublon-holon pairs if Vij > 0. This doublon-holon at¬ 
tractive correlations in the Jastrow factor play a crucial 
role for the Mott transitio n^^d^ . 


B. Correlation Factors 

By operating the correlation factors on the pair- 
product wave function, we can include many-body corre¬ 
lation effects beyond the mean-field level. In this study, 
we use the Gutzwiller factor Pg^ and the Jastrow factor 
Pj^, i.e., P = PgPj. 

The Gutzwiller factor which was introduced by 
Gutzwiller— has the form of 

Ns \ 

-g'^ni-fUiiU ( 10 ) 

where riia = cJ^Cio- and g is the variational parameter. 
The Gutzwiller factor punishes the double occupation of 
electrons on the same site in real space configurations. 


C. Quantum-Number Projections 

In general ground states of finite quantum systems, the 
symmetries of the Hamiltonian must be preserved even 
if the symmetry-breaking occurs in the thermodynamic 
limit. Furthermore, such symmetries must be preserved 
even after time evolutions as long as external fields do 
not break the original symmetries. However, conven¬ 
tional trial wave functions often break symmetries of the 
Hamiltonian. 

The quantum-number projection enables us to recover 
the symmetries of the trial wave function^. Here, we in¬ 
troduce two quantum-number projections: the spin pro¬ 
jection £‘® and the momentum projection . The spin 
projection £‘® is the projection onto the state with a total 
spin S and z-component of spin = 0. This projection 
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has a form of the integration over spin space: 

2S +1 f 

J d^]P5(cos/3)7^(^]), (13) 

where fl = (a, /3, 7 ) is the Euler angle, P 5 (cos /3) is the 5'- 
th Legendre polynomial, and TZ{V,) = is 

the rotational operator. When the one-body part \cj)) and 
the real space configurations {cc} satisfy the condition of 
= 0 , we can omit the integrations over 7 and a as 
follows: 

!</') = J df^Ps(cos/3)e“^%'^^“e*'>'^^ |<(.) 

= / d/3 sin/3Ps(cos/3) \4>). {lA) 

X 

We can efficiently estimate the integration over /3 by the 
Gauss-Legendre quadrature^. The momentum projec¬ 
tion is the projection onto the state with a total 
momentum K. This projection has a form: 

= (15) 

® R 

where Tr is the translation operator with the translation 
vector R. 


D. Backflow Correlations 

One efficient way to include correlation effects 
in a trial wave function is to include backflow 
correlation o^^d^i^^~ —. Recently, Tocchio et al. proposed 
a way of introducing backflow correlations into a Slater 
determinant for lattice models and found that the back- 
flow correlations substantially improve ground-state en¬ 
ergy of frustrated electronic systems in the region above 
intermediate strength of coupling^SiiE. In a way similar 
to the approach by Tocchio et al., the backflow correla¬ 
tions for lattice models can be implemented in the pair- 
product wave functions with the momentum projection 
as 

I/) = 

X 

X ^ R 

where \4>^) and FfX^{x -r) represent the pair-product 
wave function with backflow correlations and the Pfaffian 
of the skew-symmetric matrix X^{x-r), respectively. 
X^{xr) is defined as 

^nmi^R.) = ~ (^^ h )- ( 17 ) 

Here, xr represents the real space configuration which 
is created by shifting a configuration a: by a translation 
vector R. The site of the n-th (m-th) electron is rep¬ 
resented by in (im)- TR{in) is the site characterized by 


the position vector -|- R. For simplicity, we do not 
consider the momentum projection in the following equa¬ 
tions. The pairing amplitude with backflow correlations 
is defined by 


t,t' 


where {??(()(/} represent variational parameters. 
taken over r (r') that satisfies the following condition : 
0 < |(5| < where (5 = {S = ri^+r'-rij). 

We usually choose as the range of hopping terms in 
the Hamiltonian. Here, we drop the electron indices n 
and m for simplicity of notation, and we define {x) 

as 


^ti+r{x) = (19) 

0i^+,((r) = , (20) 

0 j ( 3 ::) — CT^z-t-r, —crhi-j-rcr)^, 5 ( 21 ) 

0z^z+7-(:r) — (I3zUz-|_z-^_crhz-|_z-cr -j- Riahij — a Ri+r') ;(22) 

where (jqi+r represents the Kronecker delta, hia- = l — riia- 
and (H)^ = (a:|H|a:) / {x\x). We impose that has the 
inversion symmetry and is independent of spin indices, 
namely Ty)))), = ryO")) Furthermore, ry™, 

is replaced with 1 when 0 )'^_|_^(a:) = 0 for any r and a. 
Note that the introduction of backflow correlations make 
computational costs heavy because we need to recalcu¬ 
late the element of pairing wave function whenever we 
generate a candidate of the next sample and calculate 
expectation values of off-diagonal operators. 

Finally, we explain the difficulties in operating the spin 
projection on the pair-product wave function with back- 
flow correlations. From the definition of the pairing am¬ 
plitudes with backflow correlations, the trial wave func¬ 
tion with the spin projection and backflow correlations is 
described as 

I/) = E / c?/3sin/3Ps(cos/3) |/) 

X 

oc E |2^) J dPsmf3Ps{cosP){x\e^^^''\x')PiX'^{x’). (23) 

x,x' 

Here, we need to take the summation over real space 
configurations x' because the skey-symmetric matrix de¬ 
pends on a:'. However, we cannot usually take this sum¬ 
mation because the number of x' grows exponentially as 
the system size increases. Thus, it is difficult to operate 
the spin projection on a trial wave function with backflow 
correlations. 


IV. RESULTS 

In this section, we show the accuracy and efficiency 
of the t-VMC method. For benchmark tests, we con- 
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sider the quench dynamics in fermionic Hubbard mod¬ 
els. We compare the t-VMC results with numerically 
exact results obtained by calculating the formal solution 
of the time-dependent Schrbdinger equation and time- 
dependent DMRG. 


A. Model and Setting 


For benchmark tests, we consider the fermionic Hub¬ 
bard model which is defined as 

T~L{t) = ~thop ^ + U{t) 'y ^ nj-friii, (24) 

<2J>,(7 i 


where thop represents the hopping amplitude between 
nearest-neighbor sites and U{t) represents the time- 
dependent onsite interaction. In this paper, we set 
^hop = 1-0. In two-dimensional cases, we treat the square 
lattice. We consider the linear-ramp quench where the 
strength of interaction U{t) is linearly changed from Ui 
to Uf during time tgi 


Uit) = 


Ui + (0 < t < tq) 

Uf {t > tq), 


(25) 


where Ui and Uf represent the strength of interaction 
before and after the quench protocol, respectively. 

For equilibrium systems, there exist many theoretical 
studies on the Hubbard models^i^^— . In the present 
models, the ground state at half-filling is believed to be 
the antiferromagnetic insulator for any nonzero 17/thop 
due to the Fermi surface nesting. Away from half-filling, 
this model shows rich phases such as correlated met¬ 
als, antiferromagnetic metal and d-wave superconducting 
state o^^i^^~ — . According to the VMC calculations, such 
superconducting states appear above C//thop ^ 6.0 in an 
interval of the doping concentration^. 

In all the t-VMC calculations, we impose a boundary 
condition so that the closed shell condition is satisfied 
for U{t) = 0.0. We choose the discrete time step At = 
1.0 X 10“^/[/(t) for U{t) > thop and At = 1.0 x 10“^/thop 
for U (t) < thop in the Runge-Kutta procedure. 


B. Comparison with Exact Results 

In this section, we mainly show the t-VMC results com¬ 
pared with “exact” results. Here, the “exact” results 
mean those obtained by directly calculating the formal 
solution of the time-dependent Schrodinger equation: 

= E (^)”w^)) + ^(^^n( 26 ) 

where the wave function is completely expanded in the 
full Hilbert space. The errors in this calculation arise 
only from the time discretization with the amplitude of 


0{At^). To preserve the total energy and the norm of 
the wave function with good accuracy, we have to choose 
a small time-grid At and a high order M. In this study, 
we choose At = 1.0 x 10“^ and M = 3 for systems at 
half filling and At = 1.0 x 10“^ and M = 5 for doped 
systems. Although the exact results can be obtained ir¬ 
respective of the interaction strength, the system size is 
severely restricted to small systems because of the expo¬ 
nentially growing computational cost with the increase 
of the system size. In order to check whether our trial 
wave function is accurate even for larger systems, we here 
compare the t-VMC results with time-dependent DMRC 
results which are expected to be highly accurate in the 
thermodynamic limit. 

In this section, we mainly show the results ob¬ 
tained by using a trial wave function with all quantum- 
number projections but without backflow correlations 
IV^(t)) = \4>it))- Although we also tried 

a trial wave function with backflow correlations |'0(t)) = 
l^^(t)) in our benchmarks, we did not observe 
any clear improvement at least for small systems. Figure 
[1] shows an example of the time evolution of the dou¬ 
ble occupancy for the one-dimensional Hubbard model 
at half Hlling and {Ui,Uf,tq,Ns)={0.0, 12.0, 5.0, 10). 
This result indicates no effect of the backflow correla¬ 
tions. Note that backflow correlations may improve the 
description of dynamics in large systems with geometri¬ 
cal frustrations in a strong coupling region because these 
correlations are important for the description of ground 
states for such system a^^d^ . 


1. One- and two-dimensional Hubbard model at half-filling 


In this section, we calculate the time evolution of sev¬ 
eral physical quantities. The physical quantities we have 
measured are the double occupancy d{t), the momentum 
distribution n{t-, k) and the spin structure factor S's(t; q) 
defined as 


Ns 


d{t) = jfYl ’ 

^ i 


(27) 


Ns 






respectively. Here, k and q are wavenumbers in the 
Brillouin zone. In addition. Si = 1/2 o-'‘'L'^o'.CT'Cicr', 

where cr are the Pauli matrices. 

Figures [5] and [3] show the time evolution of several 
quantities in one dimension for {Ui,tq, Ns)={0.0, 5.0, 16) 
and in two dimensions for {Ui,tq, Ns)={0.0, 5.0, 4x4), 
respectively. The dimension of the Hilbert space in these 
systems is NsC-n XNs Cn ~ 10®. However, our trial wave 
function has only several hundreds parameters. Never¬ 
theless, the t-VMC results well reproduce the exact re¬ 
sults. These results show that our trial wave function 
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FIG. 1. (Color online) Time evolution of double occu¬ 
pancy d{t) for linear-ramp quench in one-dimensional Hub¬ 
bard model at half-filling. The parameter set is chosen as 
{Ui,Uf,tq, Ns) = {0.0, 12.0, 5.0, 10). Red open squares rep¬ 
resent the results computed by using the trial wave function 
with the momentum projection but without backflow correla¬ 
tions and blue solid squares represent the results obtained by 
using the trial wave function with backflow correlations and 
the momentum projection. Black curve represents the exact 
result. The range of backflow parameters is chosen as 

1. To verify the effect of backflow correlations, we here opti¬ 
mized trial wave functions directly without the Monte Carlo 
integration. 


offers a highly accurate and efficient description of quan¬ 
tum dynamics for strongly correlated electron systems. 

Here, we briefly comment on a dynamical transition in 
the Hubbard models at half-filling. Several works showed 
that the dynamics of the jump An(t) is different de¬ 
pending on the strength of interaction^^— . For weak 
interactions, the jump An(t) decreases gradually from 
An(t) = 1.0 to a constant. On the other hand, for strong 
interactions, the jump An{t) exhibits a collapse-and- 
revival oscillation. In Figl^Jc), we observe clear collapse- 
and-revival oscillations in one-dimensional system espe¬ 
cially at Uf = 8.0. However, FiglSKc) shows that only 
weak oscillations are detected even at large Uf in the two- 
dimensional system. The difference appears to show a 
qualitative difference between one- and two-dimensional 
systems. However, it is difficult to conclude whether a dy¬ 
namical transition to the collapse-and-revival oscillation 
happens since the system sizes are too small to measure 
physical properties right at the Fermi surface, especially 
in two-dimensional case. To study nonequilibrium prop¬ 
erties such as the dynamical transition, we need to treat 
larger system size. We leave its analysis for a future 
study. 

Next we show the dependence on time-dependent trial 
wave functions. In Fig. IDJa), we present the t-VMC re- 



FIG. 2. (Color online) Time evolution of (a) energy per site 
E/Ns{t), (b) double occupancy d{t), (c) jump of momentum 
distribution at Fermi energy An{t) = n{t\ q = '7r/2) — n{t\ q = 
7r/2-|-7r/Ns), and (d) spin structure factor Ss(t; q — tv) for the 
linear-ramp quenches in one-dimensional Hubbard model at 
half filling. The parameter set is chosen as (Hi, tq, Ns)=(0.0, 
5.0, 16). Symbols and curves represent the t-VMC results 
and the exact results, respectively. Error bars indicate the 
statistical errors arising from the Monte Carlo sampling, but 
most of them are much smaller than the symbol sizes here 
and in the following figures. 



FIG. 3. (Color online) Time evolution of (a) energy per 
site E/Ns{t), (b) double occupancy d{t), (c) jump of mo¬ 
mentum distribution at Fermi energy An{t) = n{t\ (vr,0)) — 
n(t; (ti,h/Ns)), and (d) spin structure factor Sa{t',q = (7r,7r)) 
for the linear-ramp quenches in Hubbard model on square 
lattice at half filling. The parameter set is chosen as 
(Hi, tg, Ns)=(0.0, 5.0, 4x4). Symbols and curves represent 
the t-VMC results and the exact results, respectively. 
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suits for several different time-dependent trial wave func¬ 
tions for the quench from Ui = 0.0 to Uf = 4.0. We note 
that \(j)) instead of \(j){t)) represents the one-body part 
fixed through the time evolution at = 0)) optimized 
in the ground state before quenching. In two of the trial 
wave functions, we optimized all variational parameters 
with the momentum projection and the results agree with 
the exact ones. Other results are obtained without op¬ 
timizing some part of variational parameters or without 
operating the momentum projection. As seen in this fig¬ 
ure, these results do not reproduce the exact results with 
substantial discrepancies. Especially, we find that the 
result obtained by optimizing only the Gutzwiller fac¬ 
tor clearly disagrees with the exact one. The main rea¬ 
son for this disagreement is that Gutzwiller-type wave 
function cannot describe insulating states because of a 
lack of long-range off-site correlations which the Jastrow 
factor includes. In fact, the result obtained by optimiz¬ 
ing only the Gutzwiller-Jastrow factor shows qualitative 
agreement with the exact one. This tendency is similar 
to that in equilibrium systems where a Gutzwiller-type 
wave function fails in reproducing the physical proper¬ 
ties of the Hubbard models in finite dimensions^ldi. One 
might think that variational parameters in one-body part 
are not so important to describe the dynamics qualita¬ 
tively because optimizing the one-body part affects only 
the amplitude of the oscillation after the quench in FigH] 
(a). However, as we see clearly in Sec lIVBSl variational 
parameters in one-body part play an important role even 
in a qualitative description of nonequilibrium states. In 
Fig. Ha), we do not see any improvements by oper¬ 
ating the spin projection to \(t>{t)). However, 

for the quench to strong interaction Uf = 8.0 in Fig. 
Hb), the difference between the two wave functions with 
the momentum projection is more obvious than that for 
Uf = 4.0. These results suggest that it is better to op¬ 
erate both of the quantum-number projections on the 
trial wave function especially in the strong coupling re¬ 
gion. These results show that in order to obtain accurate 
results by the t-VMC method, we should operate the 
quantum-number projection and optimize all the varia¬ 
tional parameters. This accuracy sensitive to trial wave 
functions is similar to that for ground states (see Ref. [ni¬ 
ls^ and AnnendixlBll. 


2. Size dependence 


In the previous subsection, we have shown the bench¬ 
mark results for small systems at half-filling. In or¬ 
der to check the accuracy of our trial wave function 
for larger systems, we here compare the t-VMC results 
with highly accurate DMRG result in the thermodynamic 
limit. Since the VMC method offers the results for finite- 
size systems, we need to investigate the system-size de¬ 
pendence of our results. 

In Fig. H we present the t-VMC results of An{t) 
for sudden quench protocol (tq = 0.0) from Ui = 0.0 



FIG. 4. (Color online) Time evolution of the jump of the mo¬ 
mentum distribution near Fermi energy An{t) for the linear- 
ramp quenches to (a)t// = 4.0 and {h)Uf = 8.0 in one¬ 
dimensional Hubbard model at half-filling. The parameter 
set is chosen as {Ui,tq, Ns) = {0.0, 5.0, 16). Symbols repre¬ 
sent the t-VMC results obtained by using different trial wave 
functions. In the legend, time-dependent part of the trial 
wave functions include variational parameters we optimized. 
Solid curves represent the exact results. 


to Uf = 1.0 in the one-dimensional Hubbard model at 
half-filling. In the t-VMC calculations, we operated only 
the momentum projection to reduce numerical cost since 
spin projection is not so important in the region of small 
interaction (See Fig. H- To check the system-size de¬ 
pendence, we showed three t-VMC results with different 
system sizes. For comparisons, we also show the results 
obtained by DMRG, DMFT and DCA— . As an impurity 
solver of the DMFT and DCA calculations, the iterative 
perturbation theory (IPT) was employed. In the DCA 
calculation, the reciprocal wavevector K satisfies the fol¬ 
lowing condition: K = 2mr/Nc, where n represents inte¬ 
ger and Nc represents a cluster size. Here, Vc = 64 . As 
seen in this figure, the jump An{t) obtained by DMRG 
relaxes slowly with an oscillation. This feature is ob¬ 
served in both the DCA and t-VMC results although the 
DMFT result does not show the oscillation clearly after 
t > 1.5. However, the DCA result shows a clear deviation 
from the DMRG result at long time even when large clus- 












ter size is used. Tsuji and his coworkers have reported 
that this deviation comes from the quantum corrections 
from higher-order diagrams neglected in the IPT— . On 
the other hand, our t-VMC results have a strong size- 
dependence but approach the result of DMRG at long 
time as the system size increases. These results imply 
that our trial wave function in the t-VMC method have 
successfully included the quantum fluctuation beyond the 
DC A result. 


^n{i) 


0.86 


VMCNs=24 -Q- DMFTNs=oo 
Ns=40 —2^ DCA Ns=oo- 
Ns=48 -o- DMRGNs=oo- 
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FIG. 5. (Color online) Time evolution of jump of momen¬ 
tum distribution near Fermi energy for sudden quench from 
Ui = 0.0 to Uf = 1.0 by using several methods. The DMRG, 
DMFT and DCA results are taken from Ref. [2^. Symbols, 
dotted line, dashed line and solid curve represent the results 
obtained by t-VMC, DMFT, DCA and DMRG, respectively. 
The cluster-size employed in the DCA calculation is Ac = 64. 


3. Superconducting correlation in hole-doped Hubbard 
model in two dimensions 

The doped Hubbard model on the square lattice is one 
of the simplest models for studying the high-Tc super¬ 
conductivity in copper oxides. In this model for equilib¬ 
rium, many theoretical studies have proposed the exis¬ 
tence of d-wave superconducting states in the hole-doped 
region with strong on-site interactio n^^i^^~ —. Recently, 
some experiments have reported that nonequilibrium su¬ 
perconducting states in copper oxides have been realized 
even at room temperatur o^^d^ . In order to theoretically 
determine whether time-evolved states realize supercon¬ 
ducting state, it is necessary to calculate time evolutions 
of pairing correlation functions. 

As a benchmark for the da; 2 _j, 2 -wave pairing correla¬ 
tion functions Pd{t\r), we consider the hole-doped Hub¬ 
bard models in two dimensions. The pairing correlation 


function is defined as 

N, 


Pd{t\r) = ^ \^{Al{r,)Ad{r, + r)) 

^ f'i 

+ {Ad{ri)A\{r,+r)) . (28) 

Here, Ad{ri) represents the d^^-yi-wave superconducting 
order parameter dehned as 

^d{ri) = ~ - CiiGt)’ (29) 

^ j 

where 


Idi ^) — ~ Sr^fi{Sry,l + (30) 

is the form factor which describes the dx'^_y 2 -wave sym¬ 
metry and r = {r^,ry). 

In Figs. inKa)-(d), we compare the t-VMC results with 
the exact results of max\Pd{t; r)\ at four different t’s in 
the linear-ramp quench from Ui = 0.0 with tq = 10.0. 
Here, max|P(i(t; r)| denotes the maximum absolute value 
of pairing correlation functions \Pd(t; r)| among the same 
r = |r|. As shown in Figs. [nKa)-(d), our t-VMC results 
show good agreements with the exact results for all the 
distances at each time. This accuracy of the supercon¬ 
ducting correlations is the same as that for the ground 
state (see Appendix IBI). 

Figure (HKe) shows the dependence of max\Pd{t; r)\ on 
trial wave functions at long time t = 50.0 for Uf = 8.0. 
We again note that \(j)) instead of |(/i'(t)) represents the 
one-body part fixed through the time evolution at \(j){t = 
0)) optimized in the ground state before quenching. As 
seen in Fig. me), only the correlation function at the 
largest distance obtained by using PG{t)Pj{t) \ (j)) shows 
a large deviation from the exact result, i.e., its value is 
one order of magnitude lower than the other ones. This 
result implies that only optimizing the correlation factors 
is insufficient in reproducing pairing correlations in time 
evolution. In fact, by optimizing the amplitudes of sin¬ 
glet pairings fij in one-body part in time evolution, the 
t-VMC result at the largest distance well reproduces the 
exact one. It is important to obtain the long-range part of 
the u].ax\Pd{t; r)\ accurately because it enables us to de¬ 
tect the emergence of the superconducting phase in large 
systems. Therefore, to describe different nonequilibrium 
states flexibly, it is crucial to optimize the one-body part. 

From these benchmark results, even the superconduct¬ 
ing correlation can be well reproduced by using our trial 
wave function, which may inspire studies along this line 
in the t-VMC method. One of the intriguing studies is on 
the influence of strong laser pulse on superconductivity 
in correlated electron systems. Some recent works have 
shown that the hopping amplitude is reduced by applying 
strong laser and the relative strength of interaction to the 
transfer become effectively larger than beforeSi^— . By 
using this effect, the d-wave superconductivity may grow 
because, for ground states, it grows as the on-site inter¬ 
action increases. However, in order to confirm whether 
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nonequilibrium states show a true long-range order or 
not, calculations for larger systems are required. Stud¬ 
ies on nonequilibrium superconducting states in larger 
systems will be reported elsewhere. 


max|P^(^;r)| 




FIG. 6. (Color online) Time evolution of superconducting cor¬ 
relation functions max|Pd(f;r)| for the linear-ramp quenches 
in two-dimensional doped Hubbard model on square lattice. 
The parameter set is chosen as Ns, N) = {0.0, 10.0, 4x4, 

12). The results are measured at time (a) t — 1.0 (after 
quench protocol starts) , (b) t = 5.0 (in the middle of quench 
protocol), (c) t = 10.0 (after quench protocol ends) and (d) 
t = 50.0 (in long-time limit). The dependence on trial wave 
functions is also shown in (e) for Uf = 8.0 at f = 50.0. The 
lattice spacing is used as a unit of distance. Open symbols 
represent the t-VMC results and the other symbols represent 
the exact results. The accurate results are obtained only when 
we use lipit)) = 


V. SUMMARY AND OUTLOOK 

In summary, we have developed a time-dependent vari¬ 
ational Monte Carlo method (t-VMC) for strongly cor¬ 
related electron systems and have shown the benchmark 
results for the fermionic Hubbard model out of equilib¬ 
rium. By comparing our t-VMC results with the exact 


results, we found that our trial wave function well repro¬ 
duces exact time evolutions in both one and two dimen¬ 
sions. These results show that our trial wave function 
offers an accurate and efficient description of nonequilib¬ 
rium states in strongly correlated electron systems. 

One of the advantages of the VMC method is its wide 
applicability. In fact, the VMC method can be applied 
to complicated ab initio effective models derived by the 
downfolding scheme and contributed to identifying mech¬ 
anism of physical properties in real materials^i^. Appli¬ 
cations of the t-VMC method to such models are intrigu¬ 
ing future issues. Furthermore, the VMC method can be 
applied not only to purely electronic systems but also to 
electron-phonon coupled systems^. Therefore, it would 
be possible to study the phenomena of relaxation pro¬ 
cess and photoinduced phase transitions through phonon 
modes in real materials such as copper oxides. 


ACKNOWLEDGMENTS 

The code was developed based on the VMC code im¬ 
plemented for electron systems with contributions by D. 
Tahara and S. Morita. The authors thank T. Misawa 
for useful comments on the present study and for pro¬ 
viding them with the code for solving time-dependent 
Schrodinger equation. In the code for solving time- 
dependent Schrodinger equation, some routines in TIT- 
PACK version 2 was partly used. K.I. also thanks Y. 
Yamaji and K. Takai for fruitful discussions. The au¬ 
thors thank the Supercomputer Center, the Institute for 
Solid State Physics, the University of Tokyo for the facil¬ 
ities. This work was financially supported by Japan So¬ 
ciety for the Promotion of Science through Program for 
Leading Graduate Schools (MERIT), the MEXT HPCI 
Strategic Programs for Innovative Research (SPIRE), 
RIKEN Advanced Institute for Computational Science 
(AICS) through the HPCI System Research Project (un¬ 
der Grants No. hpI30007, hpI402I5 and hpI502II) and 
the Computational Materials Science Initiative (CMSI). 
This work was also supported by a Grant-in-Aid for 
Scientific Research (No. 22I040I0 and No. 22340090) 
from Ministry of Education, Culture, Sports, Science and 
Technology, Japan. 


Appendix A: Time-Dependent Variational Principle 
for Norm-Independent Dynamics 

In this appendix, we review the TDVP for norm- 
independent dynamics in terms of the principle of least 
actio n^^d° . We apply the principle of least action to an 
action S = J dtL{a, a) on the manifold A4 of a trial wave 
function |'(/'a)- Here, the Lagrangian L{a, ct) is described 
as 

L{a,a) = ^ ((■felV'a) - (V^IV'a)) - , 
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where ot and a represents variational parameters and 
its complex conjugates, respectively. Although the norm 
of the wave function is preserved under the exact time 
evolution, the evolution on the manifold M. may break 
the norm-conservation. To remove the restriction on the 
norm, the norm-dependent Lagrangian L(a,a.) should 
be normalized. Thus, the modified Lagrangian L{a,ck) 
for norm-independent dynamics is introduced: 


L(a, a) = L{a, a.)/ 

_ i (tWlV’g) - {'ihx\'>Pcc) _ (•fel’HlVjg) 

2 (^I'0a) 


The variation of the corresponding action 5S with respect 
to variations of parameters {'<pa \ (V'al + \ is given 

by 


5S{ol, a) 





(tpallpa) 







Stationarity of the modified action 6S = 0 leads to the 
variational equation on the manifold M 



/d-n] 

V (ihxli’a) ) 

dt 


|V'c«)=0. (Al) 


Thus, the condition of minimizing the modified action for 
norm-independent dynamics is equivalent to Eq. (ED in 
the full Hilbert space. Based on Eq.Q or Ea. dAll) . we 
can easily derive the Euler-Lagrange equation described 
as 


— 


dak 

dt 


Np 

I 


(A2) 


where a matrix S and a vector g are described as 


= aif air (V'alV'a) , gk 


daj: 


Appendix B: Dependence on Trial Wave Functions 
for Ground States of Hubbard Models 

In this appendix, to gain insight into the dependence 
of physical properties on trial wave functions in the 
nonequilibrium states, we compare them with the bench¬ 
marks for the ground states of the Hubbard models. 


Tables U and El show how physical properties depend 
on trial wave functions in the ground state of the Hub¬ 
bard models at and away from half-filling, respectively. 
To show the accuracy of our VMC results, we also show 
the results obtained by the exact diagonalization (ED) 
method. In these Tables, and \(j)opt) represent the 
Fermi sea state and optimized pair-product wave func¬ 
tion, respectively. In all the cases in Table HI there 
are large discrepancies, especially in the jump An, be¬ 
tween the results obtained by using the Gutzwiller-type 
wave function (GWF) Vg |<(>f) and those of the ED 
method. The main reason for this is that GWF can¬ 
not describe insulating states as we described in Sec IHIl 
By operating the Jastrow factor, the VMG results at 
half-filling are qualitatively consistent with those of the 
ED method. For both half-filled and hole doped models, 
the energies obtained by using our best trial wave func¬ 
tion C^^^VgVj |<^opt) agree with those of the ED 
method, and the relative errors are less than 0.5%. Figure 
[7] shows the pairing correlations obtained by using differ¬ 
ent trial wave functions in the doped Hubbard model on 
the square lattice for C//thop = 8.0, n = 12/16. As shown 
in Fig[71 j |</opt) has the best accuracy 

of the d^ 2 _y 2 -'wave superconducting correlations for all 
the distances. These trends on trial wave functions are 
similar to those in the t-VMC method. 


TABLE I. Comparison of physical quantities obtained by the 
exact diagonalization (ED) method with those by using differ¬ 
ent trial wave function for one-dimensional Hubbard model at 
half-filling. E/Ns, An and Sa{iT) represents energy per site, 
jump of the momentum distribution near the fermi energy, 
and the spin structure factor, respectively. The numbers in 
parentheses denote the statistical errors in the last digits. 



E/Ns 

An 

^^(Tr) 

t//thop=4.0, A,=16 

Vg 

-0.5280(5) 

0.843(1) 

0.490(2) 

VgVj |0f) 

-0.555(4) 

0.5257(17) 

0.677(4) 

VgVj |<^opt} 

-0.5674(4) 

0.3868(17) 

0.6878(23) 

C^^°CK-°VgVj Ifliopt) 

-0.57605(1) 

0.4220(5) 

0.7329(4) 

Exact (ED) 

-0.57660 

0.4326 

0.7277 

t//thop=8.0, A,=16 

Vg |0f) 

-0.217(2) 

0.447(3) 

0.813(5) 

VgVj |(/)f) 

-0.3170(4) 

0.168(3) 

0.904(3) 

VgVj |(/>opt) 

-0.3238(3) 

0.150(1) 

0.898(3) 

C^^°CK-°VgVj |<^opt> 

-0.32857(2) 

0.1546(3) 

0.9561(7) 

Exact (ED) 

-0.32904 

0.1578 

0.9556 
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TABLE 11. Comparison of energy per site obtained by the 
ED method with those by using different trial wave function 
for doped Hubbard model on square lattice for t//thop=8.0 at 
n=12/16. The numbers in parentheses denote the statistical 
errors in the last digit. 

E/Ns 

VgVj |</>f> -0.9373(1) 

VgVj |<(>opt) -0.9452(2) 

C^=°C'^=°VGrj |<(.opt) -0.9728(1) 

Exact(ED)-0.9774 



r 


FIG. 7. (Color online) Superconducting correlation func¬ 
tions max|Pd(r) | in the ground state of two-dimensional doped 
Hubbard model on square lattice for t//thop = 8.0, n = 12/16. 
Here, V = VgVj. The lattice spacing is used as a unit of dis¬ 
tance. Open symbols represent the VMC results and the dots 
represent the ED results. 
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